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ABSTRACT. Functionally graded materials (FGMs) are two-phase composites with continuously changing 
microstructure adapted to performance requirements. Traditionally, the overall behavior of FGMs has been 
determined using local averaging techniques or a given smooth variation of material properties. Although 
these models are computationally efficient, their validity and accuracy remain questionable, since a link with 
the underlying microstructure (including its randomness) is not clear. In this paper, we propose a numerical 
modeling strategy for the linear elastic analysis of FGMs systematically based on a realistic microstruc- 
tural model. The overall response of FGMs is addressed in the framework of stochastic Hashin-Shtrikman 
variational principles. To allow for the analysis of finite bodies, recently introduced discretization schemes 
based on the Finite Element Method and the Boundary Element Method are employed to obtain statistics of 
local fields. Representative numerical examples are presented to compare the performance and limitations 
of both schemes. To gain insight into similarities and differences between these methods and to minimize 
technicalities, the analysis is performed in the one-dimensional setting. 



1. Introduction 

Generally speaking, the ultimate goal of every design is a product which fully utilizes properties of 
materials used in its construction. This philosophy, in its vast context, naturally leads to an appearance 
of multi-phase composites with microstructure adapted to operation conditions; e.g. [Petrtyl et al., 1996, 
Bends0e and Sigmund, 2004, Ray et al., 2005]. Functionally graded materials (FGMs) present one im- 
portant man-made class of such material systems. Since their introduction in 1984 in Japan as barrier 
materials for high-temperature components, FGMs have proved to be an attractive choice for numer- 
ous applications such as wear resistant coatings, optical fibers, electrical razor blades and biomedical 
tools [Neubrand and Rodel, 1997, Uemura, 2003]. To provide a concrete example, consider a microstruc- 
ture of Al203/Y-Zr02 ceramics (Figure 1) engineered for the production of all-ceramic hip bearings. In 
this case, controlled composition and porosity allow to achieve better long-term performance and hence 
lower clinical risks when compared to traditional metallic materials [Lukas et al., 2005]. 

As typical of all composites, the analysis of functionally graded materials is complicated by the fact 
that the explicit discrete modeling of the material microstructure results in a problem which is intractable 
due to huge number of degrees of freedom and/or its intrinsic randomness. As the most straightforward 
answer to this obstacle, models with a given smoothly varying material data are often employed. When 
the spatial non-homogeneity is assumed to follow a sufficiently simple form, this premise opens the route 
to very efficient numerical schemes, such as specialized finite elements [Santare and Lambros, 2000], 
boundary element techniques [Sutradhar and Paulino, 2004], meshless methods [Ching and Chen, 2007] 
or local integral equations [Sladek et al., 2005]. Thanks to their simplicity, these methods can be rather 
easily generalized to more complex issues such as coupled thermal- mechanical problems [Noda, 1999] or 
crack propagation [Sekhar et al., 2005]. Although this approach is very appealing from the computational 
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Figure 1 . Graded microstructure of Al203/Y-Zr02 ceramics (Courtesy of J. Vleugels, 
K.U. Leuven). 



point of view, its validity remains rather questionable as it contains no direct link with the underlying 
heterogeneous microstructure. 

One possibility of establishing such a connection is to assert that the FGM locally behaves as a ho- 
mogeneous composite characterized by a given volume fraction distribution and use well-established 
local effective media theories; see, reviews by [Milton, 2002, Bohm, 2005] for more details. Local av- 
eraging techniques have acquired a considerable attention due to their simplicity comparable with the 
previous class of models; see, e.g. [Markworth et al., 1995, Cho and Ha, 2001] for an overview and com- 
parison of various local micro-mechanical models in the context of FGMs. An exemplar illustration of 
capabilities of this modeling paradigm is the work [Goupee and Vel, 2006] which provides an efficient 
algorithm for FGMs composition optimization when taking into account coupled thermo-mechanical ef- 
fects. Still, despite a substantial improvement in physical relevance of the model, local averaging meth- 
ods may lead to inaccurate results. This was demonstrated by systematic studies of [Reiter et al., 1997] 
and [Reiter and Dvorak, 1998], which clearly show that the local averaging technique needs to be adapted 
to detailed character of the microstructure in a neighborhood of the analyzed material point. When con- 
sidering the local averaging techniques, however, such information is evidently not available as all the 
microstructural data has been lumped to volume fractions only. 

Another appealing approach to FGM modeling is an adaptive discrete modeling of the structure. In 
order to avoid the fully detailed problem, a simplified model based on, e.g., local averaging techniques 
is solved first. Then, in regions where the influence of the discreteness of the microstructure is most 
pronounced, the microstructure with all details is recovered to obtain an accurate solution. Such a mod- 
eling strategy has been, e.g., adopted in [Grujicic and Zhang, 1998] when using the Voronoi cell finite 
element method introduced by [Ghosh et al., 1995] or recently in [Vemaganti and Deshmukh, 2006] in 
the framework of goal-oriented modeling. Without a doubt, this approach yields the most accurate re- 
sults for a given distribution of phases. However, its extension to include inevitable randomness of the 
microstructure seems to be an open problem. 

The systematic treatment of FGMs as random, statistically non-homogeneous composites offers, on 
the other hand, a possibility to apply the machinery of statistical continuum mechanics [Beran, 1968, 
Torquato, 2001]. In this framework, overall response of the media is interpreted using the ensemble, 
rather then spatial, averages of the involved quantities. The first class of methods stems from the de- 
scription of the material composition by a non-stationary random field. This approach was pioneered 
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by [Ferrante and Graham-Brady, 2005] and further refined in [Rahman and Chakraborty, 2007], where 
the random field description was applied to the volume fractions of the involved phases and the overall 
statistics was obtained using the local averaging methods. Such a strategy, however, inevitably leads to 
the same difficulties as in the case of deterministic analysis with a given variation of volume fractions. 
Alternative methods exploit the tools of mechanics of heterogeneous media. This gives rise to a correct 
treatment of non-local effects when combined with appropriate techniques for estimating statistics of lo- 
cal fields. Examples of FGMs-oriented studies include the work of [Buryachenko and Rammerstorfer, 2001] 
who employ the multi-particle effective field method or the study by [Luciano and Willis, 2004] based 
on the Hashin-Shtrikman energy principles; see also [Buryachenko, 2007] for a comprehensive list of 
references in this field. Both works, however, being analytically based, concentrate on deriving explicit 
constitutive equations for FGMs and therefore work with infinite bodies neglecting the finite size of the 
microstructure. 

The goal of this paper is to make the first step in formulating a numerical model which is free 
of the above discussed limitations. The microstructural description is systematically derived from a 
fully penetrable sphere model introduced by [Quintanilla and Torquato, 1997], which is briefly reviewed 
in Section 2. The statistics of local fields then follow from re-formulation of the Hashin-Shtrikman (H- 
S) variational principles introduced, e.g., in [Willis, 1977, Willis, 1981] and summarized in the cur- 
rent context in Section 3 together with the Galerkin scheme allowing to treat general bodies proposed 
by [Luciano and Willis, 2005] or [Prochazka and Sejnoha, 2003]. Section 4 covers the application of the 
Finite Element Method (FEM) following [Luciano and Willis, 2005, Luciano and Willis, 2006] and the 
Boundary Element Method (BEM) in the spirit of [Prochazka and Sejnoha, 2003]. Finally, based on re- 
sults of a parametric study executed in Section 5, the comparison of both numerical scheme when applied 
to FGMs modeling is performed in Section 6 together with a discussion of future improvements of the 
model. In order to make the presentation self-contained and to minimize technicalities, the attention is 
restricted to an one-dimensional elasticity problem (or, equivalently, to a simple laminate subject to body 
forces varying in one direction; cf. [Luciano and Willis, 2001]). 

In the following text, we adopt the matrix notation commonly used in the finite element literature. 
Hence, a, a and A denote a scalar quantity, a vector (column matrix) and a general matrix, respectively. 
Other symbols and abbreviations are introduced in the text as needed. 



2. Microstructural model 

As already indicated in the introductory part, the morphological description adopted in this work 
is the one-dimensional case of a microstructural model studied in [Quintanilla and Torquato, 1997]. A 
particular realization can be depicted as a collection on N rods of length I distributed within a structure 
of length L, see Figure 2. The position of the i-th rod is specified by the x coordinate of its reference 
point Xi, which in our case coincides with the midpoint of a rod. 



linn i in mi II 1 1 1 







Figure 2. Example of microstructural model realization. 



The microstructure gradation is prescribed by an intensity function p(x), with the product p(x) dx 
giving the expected number of reference points in an infinitesimal neighborhood around x. Using the 
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theory of general Poisson processes, the probability of finding exactly m points located in a finite-sized 
interval / is given by [Quintanilla and Torquato, 1997] 



Fm(/) = ^^exp (-Mp(J)) with ^ p (I) = Jp( 



x) dx. (1) 

Further, to provide a suitable framework for the description of microstructure related to the model, we at- 
tach a symbol a to a particular microstructure realization (e.g., Figure 2) from a sample space § endowed 
with a probability measure p. Then, the ensemble average of a random function f(x, a) is defined as 1 



</>(*) = / f(x,a)p(a)da. (2) 

Now, interpret Figure 2 as a distribution of "white" and "black" phases. For a given configuration a, 
the distribution of a phase r is described by the characteristic function Xr{x; a) 

, . f 1 if x is located in phase r, 
Xr(x;a) = { Q Qtherwise) (3) 

where r = 1 is reserved for the white phase (matrix) while r = 2 denotes the black phase (rod). The 
elementary statistical characterization of the model is provided by the one-point probability function S r 

S r (x) = ( Xr )(x) (4) 

giving the probability of finding a point x included in the phase r. Recognizing that the probability of 
locating x in the white phase coincides with the probability that the interval 

I(x) = [x-£/2,x + £/2] (5) 

will not be occupied by any reference point and using Equation (1), we obtain 

rx+e/2 \ 

p(t) dt . (6) 

/2 J 

The one-point probability function S^ix) follows from the identity 

S 1 {x) + S 2 {x) = l, (7) 



Si(aO = P Q {I(x)) = exp (- j* 



which is a direct consequence of the adopted definition of the characteristic function; recall Equation (3). 
By analogy, we can introduce the two-point probability function S rs 



S rs (x,y) = / Xr(x,a)xs(y,a)p(a)da, (8) 
Js 

quantifying the probability that a point x will be located in phase r while y stays in the phase s. For 
r = s = 1, the descriptor coincides with the probability that the union of intervals I(x) and I(y) will not 
be occupied by a reference point, yielding 

S n (x,y) = P (I(x)Ul(y)). (9) 

The remaining functions S rs can be directly expressed from Sn using relations summarized in Appen- 
dix A. 



'To simplify the exposition, we introduce the following notation: for a real-valued random function f(x, a) : R x § — » R, 
by writing f(x; a) we mean a deterministic function of x G R related to a given fixed realization (i.e., f(x; a) : R — > E). In 
other words, it holds f(x; a) := f(x, 0)\ 
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FIGURE 3. Examples of one- and two-point probability functions for a = 0.25L, b 
0.75L, L = 1 m, p a = -log(0.25/£) and p b = - log(0.75/£). (b) £ = 0.1. 



Finally, to provide a concrete example, consider a piecewise linear intensity function 

p a < x < a 

Pa + k p (x — a) a < x < b 

Pb b < x < L ' 

otherwise 



p(x) 



(10) 



where k p = (p — p a )/{b — a). The corresponding one- and two-point probability functions, evaluated 
using an adaptive Simpson quadrature [Gander and Gautschi, 2000], are shown in Figure 3. Obviously, 
the shape of one-point probability function directly follows from the intensity profile (up to some bound- 
ary effects due to extension of p by zero outside of and smoothing phenomena with lengthscale i 
demonstrating the "geometrical" size effect present in the model). The two-point probability function 
then contains further details of the distribution of individual constituents. 



3. HASHIN-SHTRIKMAN VARIATIONAL PRINCIPLES 

The introduced geometrical description provides a solid basis for the formulation of a stochastic model 
of one-dimensional binary functionally graded bodies. In the sequel, we concentrate on the simplest case 
of linear elasticity with deterministic properties of single components. 
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FIGURE 4. One-dimensional elasticity problem associated with realization a. 



3.1. Problem statement. Consider a bar of unit cross-section area, represented by the interval O, = 
(0,L) with the boundary = {0,L}, fixed at dfl u , subject to a body force b(x) and tractions t at 
d^l f , see Figure 4. For a given realization a, the displacement field u(x; a) follows from the energy 
minimization problem 

u(x; a) = arg min IL(v(x);a), (11) 

ti(i)ev 
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where arg imn xe x f( x ) denotes the minimizer of / on X, V is the realization-independent set of kine- 
matically admissible displacements, v is a test displacement field and the energy functional II is defined 
as 



(12) 



H(v(x); a) = — I e(v(x))E(x; a)e(v(x)) dx — f v (x)b(x) dx — (y(x)t(x))\ 
2 Jn Jn 

with the strain field e(v(x)) = 4 s (a;) and the Young modulus E in the form 

E(x;a) = xi(x;a)E 1 +X2(x;a)E 2 , (13) 

where Ei denotes the deterministic Young modulus of the i-th phase. 

Now, given the probability distribution p(a), the ensemble average of displacement fields follows 
from the variational problem [Luciano and Willis, 2005]: 

(u)(x) = / I arg min H(v(x),a) ) p(a) da. (14) 

Js \ »(i,a)eVxS / 

In theory, the previous relation fully specifies the distribution of displacement fields. The exact specifica- 
tion of the set § is, however, very complex and the probability distribution p(a) is generally not known. 
Therefore, the solution needs to be based on partial geometrical data such as the one- and two-point 
probability functions introduced in Section 2. 



3.2. Hashin-Shtrikman decomposition. Following the seminal ideas of [Hashin and Shtrikman, 1962] 
and [Willis, 1977], the solution of the stochastic problem is sought as a superposition of two auxiliary 
problems, each characterized by constant material data E°. 
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Figure 5 . Problem decomposition; (a) deterministic reference case, (b) stochastic po- 
larization problem. 



In the first "reference" case, see Figure 5(a), the homogeneous structure is subject to the body force b 
and the boundary tractions t. The second "polarization problem", shown in Figure 5(b), corresponds to a 
homogeneous body loaded by polarization stress r arising from the stress equivalence conditions: 

a(x; a) = E(x; a)e(x; a) = E e(x; a) + t(x; a). (15) 

The unknown polarization stress now becomes a new variable to be determined as the stationary point 
of the two-field Hashin-Shtrikman-Willis functional; e.g. [Willis, 1977, Prochazka and Sejnoha, 2004] 
and [Bittnar and Sejnoha, 1996, Chapter 1.8] 

(u(x; a), t(x; a)) = arg min stat U(v(x), 0(x; a); a), (16) 
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where 8 denotes an admissible polarization stress from the realization-dependent set T(a) , arg stat xg x f(x) 
stands for a stationary point of / on X and a new energy functional U is defined as 

U{v{x),9(x;a);a) = - [ e{v(x))E°e{v{x)) dx - [ v(x)b(x) dx - (y(x)t(xj) \ 9Qt (17) 



(E(x; a) — E°) 6(x; a) dx. 



The minimization with respect to v in Equation (16) can be efficiently performed using Green's func- 
tion technique. To that end, we introduce a decomposition of the displacement field 

u(x;a)=u°(x)+u 1 (x;a), (18) 

where u° solves the reference problem, while u 1 denotes the displacement field due to a test stress 
polarization field 9. Note that the determination of u° is a standard task, which can be generally solved 
by a suitable numerical technique (cf. Sections 4. 1 and 4.2). By introducing the Green function of the 
reference problem satisfying 

E°^(x,y)+6(x-y) = (19) 
with boundary conditions (n denotes the outer normal, recall Figure 5) 

G°{x, y) = for x e dtt u , T°(x, y) = E° 9G n{x) = for x G dtt\ (20) 

we relate the u 1 component and the associated strain field e 1 to the polarization stresses 9 via, 
cf. [Luciano and Willis, 2005], 

u\x;a) = - [ dG °j ) X ' y) 9(y;a)dy = - [ A°(x, y)9(y; a) dx, (21) 
Jn dy J n 

e{u\x;a)) = - [ V) 6(y; a) dy = - [ T°(x, y)9(y; a) dx. (22) 

Jn oxoy Jq 

By exploiting the optimality properties of the minimizing displacement u(x; a) and upon exchang- 
ing the order of optimization, Equation (17) can be, after some steps described in, e.g. [Willis, 1981, 
Luciano and Willis, 2005], recast solely in terms of the polarizations: 

r(x;a) = arg stat H (9(x; a); a) (23) 

0(z;a)eT(a) 

where the "condensed" energy functional is defined as 

H(9(x;a);a) = min U(v(x), 9(x; a); a) = n°(n°(x)) + / 6(x; a)e (u°(x)) dx (24) 

- I [ 9(x;a)(E(x;a)-E°y 1 8(x;a)dx-l [ [ 9(x;a)T°(x,y)9(y;a) dxdy 
2 J 2 Jq Jq 

with n° denoting the total energy of the reference structure. 

With the Hashin-Shtrikman machinery at hand, the stochastic problem introduced by Equation (14) 
can be solved by repeating the previous arguments in the probabilistic framework. In particular, taking 
the ensemble average of Equations (18) and (21) yields 

(u)(x) = u°(x) - [ A°(x, y) {r){y) dy (25) 
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where the expectation (r) is a solution of the variational problem 

(t)(x) = / I arg stat H(9(x, a), a) I p(a) da. (26) 

Js V 9(x,a)eT(a)x§ / 

Again, due to limited knowledge of detailed statistical characterization of phase distribution, the pre- 
vious variational problem can only be solved approximately. In particular, we postulate the following 
form of polarization stresses: 

t(x, a) w xi{x, a)n(x) + X2(x, a)T 2 (x), 9(x, a) ps xi(%, a)8 1 (x) + x 2 («, a)9 2 {x), (27) 

where r r and 9 r are now the realization-independent polarization stresses related to the r-th phase. Plug- 
ging the approximation into Equation (26) leads, after some manipulations detailed in e.g. [Willis, 1981, 
Sejnoha, 2000], to the variational principle 



(n(x),r 2 (x)) 



2 r 

arg stat IJ°(u°(x)) + > / 8 r (x)S r (x)e (u°(x)) dx 
6 (e 1 (x),e 2 (x)) f^J n K " 

-V / 9 r (x)S r (x) [E r — E°) l 8 r (x)dx 

[ 8r(x)S rs (x,y)T (x,y)9 s (y)dxdy, (28) 



r=l s=l 



i.e. the "true" phase polarization stresses r r satisfy the optimality conditions (r = 1, 2) 

2 

[ 8 r {x)S r {x)(E r -E°y l T r {x)dx + y\ [ [ 9 r {x)S rs {x,y)T {x,y)r s (y)dydx 
Jq s=1 JnJn 



V(x)S r (x)e(u°(x))dx (29) 



for arbitrary 9 r 



3.3. Discretization. Two ingredients are generally needed to convert conditions (29) to the finite-dimensional 
system: (i) representation of the reference strain field and the Green function-related quantities and 
(ii) discretization of the phase polarization stresses. The first step is dealt with in detail in Section 4; 
now it suffices to consider the approximations 

e°' ho (x), A°' ho (x) mdr°' ho (x,y), (30) 

where ho denotes a parameter related to the discretization of the reference problem. 2 

Next, we reduce Equation (29) to a finite-dimensional format using the standard Galerkin procedure. 
To that end, we introduce the following discretization of the phase polarization stresses 

t t {x) » N Thl (x)d T r hoh \ 9 r (x) » N Thl (x)4 h \ (31) 

where N rfel is the matrix of (possibly discontinuous) shape functions controlled by the discretization 
parameter hi; d® hl and d T r h ° hl denote the degrees-of-freedom (DOFs) of trial and true polarization 



To be more precise, the goal is not to obtain accurate estimates of the Green function-related operators themselves, but 
rather to approximate the action of the operators; see Section 5.2 for further discussion. 
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stresses, the latter related to the discrete Green function. Introducing the approximations (31) into the 
variational statement (29) and using the arbitrariness of d r hl leads to a system of linear equations 

2 

K T r hl d T r hohl + ^ Ks° hl d T s hohl = R T r hohl (32) 

with the individual terms given by (r, s = 1, 2) 

K T r hl = [ N Thl (x) T S r (x) [E r - E°] ^ N Thl (x) dx, (33) 



n 



T^Th hi 



: / / N^(x) T S rs (x,y)T°^(x,y)N Th ^y)dxdy, (34) 
Jo, Jn 

R rh hi = f ^ h ^ x ) T S r (x)e°' ho (x)dx. (35) 
Jn 

Finally, once the approximate values of phase polarization stresses are available, the elementary sta- 
tistics of the displacement field follow from the discretized form of Equation (25) 

(u)(x) « (u) h ^(x) = u°> h °(x) -J2 ( K J^°' ho (x,y)S r (y)-N^(y)dy ) j cT r h ° h K (36) 

Note that additional information such as conditional statistics can be extracted from the polarization fields 
in post-processing steps similar to Equation (36); see [Luciano and Willis, 2005, Luciano and Willis, 2006] 
for more details. 

4. Reference problem and Green's function-related quantities 

4.1. Finite element method. The solution of the reference problem follows the standard Finite Element 
procedures, see e.g. [Bittnar and Sejnoha, 1996, Krysl, 2006]. Nevertheless, we briefly repeat the basic 
steps of the method for the sake of clarity. 3 The reference displacement u° follows from the identity 

e(v(x))E°e(u°(x))dx = J v{x)b(x) dx + (v(x)t(x)) \ dQt , (37) 

which should hold for any test function v G V. Within the conforming finite element approach, the 
unknown displacement u° and the test function v together with the associated strain field are sought in a 
finite-dimensional subspace V C V 

u°(x) w u°' ho {x) = N uho {x)d uho , v(x) « v ho {x) = N uho (x)d vho , (38) 

s(u°(x)) e(u°' ho (x)) = B uh °(x)d uh °, s(v(x)) e{v h °{x)) = B uh °(x)d vh °, (39) 

where N u/l ° is the displacement interpolation matrix and B u/l ° denotes the displacement-to-strain matrix. 
Using the discretized fields, Equation (37) reduces to the system 

K uh od uh = R uh^ (4Q) 

where 

K uh = I B uho (x) T E°B uh °(x)dx, (41) 
In 



R uho = J N uho (x) T b(x)dx+ (N uho (x)t(x) 



(42) 

en* 



■3 

Recall that for simplicity, we assume homogeneous Dirichlet boundary conditions only. The treatment of the non- 
homogeneous data can be found in [Luciano and Willis, 2005, Appendix A]. 
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Solving the system for d uh ° enables us to obtain the e°' h ° approximation using Equation (39)i. 

The discretized version of the Green function follows from Equation (37) with t = and b = 5(y — x), 
cf. Equations (40) and (42), 

G°{x,y) » G°' h °(x,y) = N uh °{x) (K uh ° Y N uh °{y) T . (43) 



The remaining Green function-related quantities can now be expressed directly from Equations (21) 
and (22), leading to 

A°(x,y) fs A°' ho (x,y) =N uho (x)(K uh °y 1 B uho (y) T , (44) 



T°(x,y) ~ T°' ho (x,y) = B uh °(x) (k u M ^'^(y) 7 . (45) 

4.2. Boundary element discretization. Following the standard Boundary Element Method procedures 
(e.g. [Bittnar and Sejnoha, 1996, Duddeck, 2002]), we start from the Betti identity written for the refer- 
ence problem: 

^(0E°u°(0 dC = (n(0e(v(0)E°u°(0 ~ «(0*°(0) \ d n( ~ [ <0K0 ^ (46) 

J Q 



n 



and apply the test displacement in the form 

v(0 = G > oo (H,x), (47) 
where G 0,o ° is the infinite body Green's function defined as the solution of 



£° ^2 +<*(s-fl=0, G°>°°(S,x) =G ^(x,0- (48) 



In the one-dimensional setting, this quantity is provided by e.g. [Luciano and Willis, 2001, Eq. (13)] 

G 0,oo {x ^ ) = _±_ ]x _^ (49) 

and the integral identity (46), written for any receives the form: 

u 0M {x) = ( G ^(x,Ot ' ho (0-T ^(x,Ou OM (0) , + [ G°>°°(x,Z)b(0<%, 
where the tractions T 0,oo (a;, £) are defined analogously to Equation (20)2: 



(50) 



r u '°°(x, = E» QjH^™(0 = [H(x - - - 1 n(0 for £ e On (51) 

and H denotes the Heaviside function. Imposing the consistency with boundary data for x —y 0+ and 
x — > L- yields the system of two linear equations 

E°u°' ho (L) — E°u°' h ° (0) — Lt°' h ° (L) = [ £fc(£)d£, (52) 

Jn 

E°v°' ho (L) - E°u°' ho (0) - Lt°' h °(0) = f(L- d£. (53) 

Jn 



MICROSTRUCTURE-BASED MODELING OF ELASTIC FGMS: ONE DIMENSIONAL CASE 



11 



Since one component of the pair (u 0,h ° ,t 0,h °) is always specified on dCl and dCl u ^ 0, the previous 
system uniquely determines the unknown boundary data (i.e. u°' h ° on dCl 1 and t°> h ° on dCl u ), needed to 
evaluate Equation (50). 4 

Making use of the identity 2E°d x G 0,co (x, £) = 1 — 2H(x — £), the associated strain field can be 
expressed as 



J0,ho 



(x) 



(dG °>°°(x,Q t0M 
dx 



(0 



+ 



<9x 



-6(0 



2£° 



i '^ (L) — t u ' ft ° (0) 



0,/io/ 



j\(a)dt+ j\(o^y (54) 



Analogously to the Finite Element treatment, the expression for the finite-body Green function starts 
from Equation (46) with b = 5{y — £) and boundary data (20). Following the specific form of Equa- 
tion (50) (and allowing for a slight inconsistency in notation), we introduce a decomposition of the Green 
function into the discretization-independent infinite-body part and the discretization-dependent boundary 
contribution: 

G°(x, y) « G°'°°(x, y) + G°' h ° (x, y), (55) 
where the boundary part, written for x £ Cl and y G Cl, assumes the form 



G°' h °(x,y) = (G ^(x,OT ' h °^,y)-T ^(x,OG°' h °(^yj) 



(56) 



with the boundary displacements G 0,h ° and tractions T°' h ° at £ E 50 due to a unit impulse at y deter- 
mined from a linear system (compare with Equations (52) and (53)) 

E°G OM (L,y)-E G ' ho (0,y)-LT ' ho (L,y) = y, (57) 
E°G°' h °(L,y) - E°G 0,h °(0,y) - LT°' h °(0,y) = L-y. (58) 

Expression for A is derived following an analogous procedure. We exploit the infinite-body-boundary 
split 

A (x,y)nA >°°(x,y) + A°> h °(x,y) (59) 
and obtain the first part directly from the definition (21) 



A u <°°0r,y) 



dG°>°°(x,y) 
dy 



2E° 



(2H(x -y)-l) 



The boundary-dependent part now follows from 

A°> h °(x,y)= (G u '°°(x,0 



0iOO , .dT°> h °(t,y) _„ i00/ _ ^dG°> h °(Z,y) 



dy 



T u '°°{x,0- 



Oy 



an(o 



with the y-sensitivities of boundary data evaluated from (57)-(58): 



E' 



dG°> h °(L,y) EO dG°> h °(0,y) T dT°' h °(L, y) 



dy 



dy 



L- 



dy 



dG°' h °(L,y) EQ dG°' h °(0,y) T dT°^(0,y) 



dy 



dy 



L 



dy 



-1. 



(60) 
(61) 

(62) 
(63) 



It can be verified that Equation (50) now provides the exact one-dimensional displacement field rather than an approximate 
one. Nevertheless, to keep the following discussion valid in the multi-dimensional setting and consistent with Section 4.1, we 
keep the index "ho" in the sequel. 
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The BEM-based approach is completed by approximating T° function. In particular, we get 

T°(x,y) « r ^(x,y)+r°> h °(x,y) (64) 

r°>°°(*,y) = dA °2 X,V) =^S( X -y) (65) 



T°> h °(x,y) 



dA°' h ° (x,y) ( (x, dT°' ho (f , y) 



<9x y dx 5y 

1 fdT°' h °(L,y) dT°' ho (0,y) 



an(0 

(66) 



2£° V dy ay 
Finally note that the previous procedure can be directly translated to multi-dimensional and/or vectorial 
cases; see [Prochazka and Sejnoha, 2003, Section 3] for more details. 

5. Numerical examples 

Before getting to the heart of the matter, we start with converting the relations (32)-(36) into the fully 
discrete format by replacing the integrals by a numerical quadrature and selecting a specific form of 
shape functions N rftl . To that end, we introduce a set of integration points {(1X2, ■ ■ ■ , Cn c } as well 
as associated integration weights {wi , W2 , ■ ■ ■ , } an d evaluate the components of the system matrix 
and right-hand side vector as 

« ^«; i N Tftl (Ci) T ^r(Ci)[^r--B°] 'N^CCO. (67) 

i=i 

Ks ohi « x;x;^^N^(o) T s rs (c J ,o)r°' h °(o,o)N r?ii (o) 
i=i 3=1 

+ [ [ N^(x) T S rs (x,y)T ^(x, y )-N^(y)dxdy, (68) 
Jn Jn 

Rrhohi ra ^N^ 1 (Ci) T ^r(Ci)e 0, ' l0 (Ci) ; (69) 
i=l 

2 

r=l i=l 

2 

- E( J ( A °' 00(x ' y) ' s '- (y)Nr/ll(y)dy )^ ' 11 ' (70) 

with the convention r ' 00 = A ' 00 = for the FEM-based approximation of the polarization problem. 
The basis functions and integration schemes employed in the sequel, based on a uniform partitioning of 
SI into iV e cells Q e of length hi = L/N e , are defined by Figure 6. In particular, the specification of the 
polarization stress in terms of Po shape functions requires 2N e DOFs (i.e. one DOF per cell and phase), 
while Pi and P_i discretizations are parametrized using 2(N e + 1) or 4N e values, respectively. 

Note that the BEM-related infinite-body contributions appearing in Equations (68) and (70) are still 
kept explicit, as they are available in the closed form and can be treated separately. In the present 
case, action of the r ' 00 operator is local (recall Equation (66)), while the quantities related to A ' 00 
are evaluated at cell nodal points and linearly interpolated to the interior of a cell to account for the 
discontinuity of the integrand. 
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Figure 6. Choice of shape functions and integration points related to the e-th cell; 
(a) piecewise-constant basis functions (Pq) and the Gauss-Legendre quadrature of or- 
der 1 (GLi), (b) piecewise-linear discontinuous basis functions (P_i) and the Gauss- 
Legendre quadrature of order 2 (GL2), (c) piecewise linear continuous basis func- 
tions (Pi) and Newton-Cotes quadrature of order 1 (NCi); o cell nodes, O degrees 
of freedom, ■ integration points. 

To summarize, the following factors significantly influence the accuracy of the discrete Hashin-Shtrikman 
scheme: 

• approximation of the Green function of the comparison body, 

• basis functions and numerical quadrature used to discretize the polarization problem, 

• Young's modulus of the reference body E°, 

• contrast of the Young moduli of individual phases (E2/E1), 

• characteristic size of microstructure with respect to the analyzed domain (£/L). 

All these aspects are studied in detail in the rest of this Section. Two representative examples of structures 
subject to a uniform body force b and homogeneous mixed and Dirichlet boundary data are considered, 
see Figure 7: 

statically determinate structure : u(0, a) = 0, t(L,a)=0, (71) 
statically indeterminate structure : u(0,a)=0, u(L,a) = 0. (72) 

In both cases, the heterogeneity distribution is quantified according to the model introduced in Section 2 
with the one- and two-point probability functions plotted in Figure 3. Moreover, taking advantage of 
the one-dimensional setting, we systematically compare the obtained numerical results against reliable 
reference values determined by extensive Monte-Carlo (MC) simulations introduced next. 

5.1. Direct simulation results. For the purpose of the following discussion, the reference values of the 
average displacement fields (u)mc(x) together with the 99.9% interval estimates [(u_)mc(^), (u+)mc(x)} 
are understood as the piecewise linear interpolants of discrete data sampled by MC procedure described 
in detail in Appendix B. In addition, the homogenized displacement field ur(x), corresponding to a 
deterministic structure with the position-dependent elastic modulus 

1 S 1 (x) 

is introduced to asses the performance of the local averaging approach. Figure 7 stores several represen- 
tative results plotted using dimensionless quantities. 

As apparent from Figure 7, the obtained statistics of overall response exhibits rather narrow confi- 
dence intervals, implying the reliability and accuracy of the MC estimates. For the statically determinate 
structure, the locally homogenized solution coincides with the ensemble average of the displacement 
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FIGURE 7. Reference Monte-Carlo solution; (a) statically determinate and (b) inde- 
terminate problems; MC results correspond to 99.9% confidence interval estimates, H 
refers to homogenized solution. 



fields, as demonstrated by the overlap of simulation results with homogenized data. The converse is 
true (with the 99.9% confidence) for the statically indeterminate case, where these two results can be 
visually distinguished from each other. The mismatch (which increases with increasing EijE\ or £/L) 
clearly demonstrates that even in the one-dimensional setting the local averaging may lead to incorrect 
values when treating non-homogeneous random media. These results are consistent with the fact that 
in the statically determinate case, the stress field a(x, a) is independent of a as follows from the the 
one-dimensional equilibrium equations d x a(x, a) + b(x) = and the deterministic value of traction 
at x = L due to boundary condition provided by Equation (71)2- In the latter case, however, the trac- 
tion value as well as stress field become configuration-dependent. Such effect does not appear in the 
classical homogenization setting, where for t/L — > the harmonic average is known to represent the 
homogenized solution exactly, cf. [Murat and Tartar, 1997]. This result naturally justifies the application 
of approaches based on higher-order statistics to FGMs, with the H-S method being the most prominent 
example. 

5.2. Effect of the Green function approximation. In order to illustrate the effect of approximate 
Green's function, we restrict our attention to the statistically determinate structure and employ the stan- 
dard piecewise linear basis functions N u?i0 to evaluate T 0,h ° function in the FEM setting using Equa- 
tion (45). Figure 8 allows us to perform the qualitative assessment of the results for different choices of 
basis functions, the integration scheme and the discretization parameter ho- 

Evidently, a suitable choice of discretization parameter ho is far from being straightforward. From all 
the possibilities presented in Figure 8(a), only the combinations h\ = ho with Po/GLi discretization of 
the polarization problem and h\ = 2ho with P_i/GL2 scheme are capable of reproducing the homog- 
enized solution, while all the remaining possibilities lead to inaccurate results often accompanied by an 
oscillatory response. On the other hand, the /io-independent BEM-based solutions show correct response 
for all discretizations of the polarization problem (and are virtually independent of the scheme used due 
to sufficiently low value of hi parameter, cf. Section 5.5). 

To shed a light on such phenomenon, consider the FEM approximation of r 0,h ° function plotted 
in Figure 9(a). The piecewise linear basis functions used to express the reference displacements imply 
the piecewise constant values of F 0,h °(x, y) approximating the exact expression -^u5(x — y), cf. Equa- 
tion (64). As pointed out by [Luciano and Willis, 2006], however, the accuracy of the HS scheme is 
governed by the correct reproduction of the action of the r°(x, y) operator rather than the local values. 
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Figure 8. Influence of approximate Green's function for the statically determinate 
problem; (a) FEM-based solution, (b) BEM-based solution; E 2 /Ei = 10, E°/E 1 = 5, 
l/L = 0.1, h\/l = 0.25. 




FIGURE 9. (a) Finite element approximation of the Green's function, (b) convergence 
rates of FEM vs. BEM; E 2 jEx = 10, EP/E-l = 5, l/L = 0.1. 



In the present context, it follows from Equation (69) that such requirement is equivalent to the accu- 
rate representation of the V (x,y) operator action for x coinciding with the integration points related 
to the selected numerical quadrature. It can be verified that this condition is satisfied only for the two 
aforementioned discretizations of the reference problem. In particular, for the Fo/GL\ combination we 
obtain (see Figures 9(a) and 10(a)) 

r 1 (fhohi 

r°(c e) £K(0d£ « ^r ^(Ce,Ce)^ ofel = = (74) 



i.e. the numerical scheme reproduces the action of T exactly. 
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Using Figures 9(a) and 10(b), we find the analysis of the P_i/GL2 discretization completely analo- 
gous: 



h l _J_jrh h 1 _ a 2e-l,; 



2e-l,r 



rhohi 
l,r 



Thohi 



j^T (C2e,OTr(0<% ~ ™2eT°' h ° (C 2 e, C 2e )<f 1 = , (75) 

which explains the good performance of the particular discretization scheme. 

To allow for the quantitative comparison, we exploit the fact that the exact solution is available for the 
statically determinate case and introduce a relative L 2 error measure 

H 0hl _ ||(n)^(x)-n H (x)|| L2(0) 

'/H — II / Ml • v'OJ 

H«H(aOIU a (n) 

The resulting convergence rates of the FEM- and BEM-based approach are shown in Figure 10(b) with 
integrals in Equation (76) evaluated using an adaptive Simpson quadrature [Gander and Gautschi, 2000] 
with the relative accuracy of 10 -6 . Clearly, the performance of the BEM-based scheme is slightly su- 
perior to the (properly "tuned") FEM approach. By a sufficient resolution of the reference problem, 
however, both approaches become comparable. Moreover, the results confirm good performance of Po 
and Pi schemes when compared to the P_i discretization, which requires about twice the number of 
DOFs of former schemes for the same cell dimensions hi (recall Figure 6). Similar conclusions can also 
be drawn for the statically indeterminate case. Therefore, in view of the above comments, we concentrate 
on the BEM approach in the sequel and limit the choice of basis functions to Po and Pi only. 

5.3. Influence of the integration scheme and basis functions. Thus far, we have investigated the com- 
bination of the "polarization" numerical quadratures and shape functions, for which the location of in- 
tegration points coincides with the position of DOFs. Figure 1 1 shows the convergence plots for the 
relevant basis function/integration scheme pairs. To address also the statically determinate case, the 
relative error is now related to MC data, leading to the definition 

h 0hl _ \\(u) h0h Hx)-(u)Mc(x)\\ L2{n) 
VMC ~ \\(u)Mc(x)\\ L2m ■ (7?) 

In addition, two comparative values are introduced: the relative error of the homogenized solution H (de- 
termined by Equation (77) with (u) h ° hl replaced by «h) and the relative error associated with (u_)mc 
or (u + )mc function, appearing as the Interval Estimate (IE) line. 




Figure 10. Valid combinations of discretized Green's function and polarization 
stresses; (a) Pq/GLi, (b) P_i/GL2; □ finite element nodes. 
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Figure 1 1 . Influence of the choice of numerical discretization; (a) statically determi- 
nate and (b) indeterminate structures; E^/Ex = 5, £/L = 0.1, E°/Ex = 3; IE denotes 
the error associated with the 99.9% confidence interval estimate. 



For the statistically determinate structure, the observed behavior is rather similar to the one reported 
in Section 5.2. In particular, Figure 11(a) confirms that the H-S solution quickly reaches the accuracy 
comparable with the confidence intervals (indicated by the grey area) and eventually converges to the 
homogenized solution, with the exception of Fx/GLx combination resulting in a singular system ma- 
trix (32). Moreover, the superiority of the GL2 quadrature over lower-order scheme is evident; the 
proper representation of spatial statistics seems to be more important than smoothness of the polarization 
shape functions. 

Figure 1 1(b) shows the results for the statically indeterminate case. With 99.9% confidence, the results 
quantitatively demonstrate that the homogenized solution differs from the MC data. The H-S solution 
gives the error about 50% of the value of the homogenized solution, but ceases to attain the accuracy 
set by the confidence interval. It should be kept in mind that the H-S result actually delivers an estimate 
pertinent to the fixed value of parameter E° and all random one-dimensional media characterized by the 
two-point statistics (9). 

5.4. Influence of the reference media and phase contrast. Having identified the intrinsic limitation 
of the H-S approach, we proceed with the last free parameter of the method: the choice of the reference 
medium. To that end, we introduce the following parameterization of the Young modulus 

E° = (1 - u)E x +ujE 2 . (78) 

Note that for the phases indexed such that Ex < E2, lo = and lo = 1 correspond to the rigorous lower 
and upper bounds on the ensemble average of the energy stored in the structure and, consequently, to the 
positive- or negative-definite system matrix [Prochazka and Sejnoha, 2004, Luciano and Willis, 2005]. 
The intermediate values lead to energetic variational estimates and to a symmetric indefinite system 
matrix. Figure 12 illuminates the effect of lo, plotted for two representative contrasts of phase moduli 
and hx/i ratios. 

In the first case, see Figure 12(a), the choice of the reference media has almost negligible effect on 
the H-S solution error; the slight influence observed for the coarse discretization completely disappears 
upon cell refinement. This is not very surprising as the homogenized solution depends on the first- 
order statistics only, recall Equation (73), and as such can be retained by the discrete H-S method (up 
to controllable errors) for any choice of E°. Results for the statically indeterminate structure, on the 
other hand, show a significant sensitivity to the value of u. By a proper adjustment of the reference 
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Figure 12. Influence of the choice of the reference media; (a) statically determinate 
and (b) indeterminate structures; £/L = 0.1, P0/GL2 discretization. 



medium, the error can be reduced by an order of magnitude and eventually reach the accuracy of extensive 
MC sampling. With increasing phase moduli contrast, however, the range of such uj values rapidly 
decreases; for E2/ E\ = 100 one needs to satisfy 9 • 10 -4 < u) < 1.5 ■ 10 -3 in order to recover the MC 
results. It is noteworthy that these values agree rather well with the particular choice of reference media 
used by [Matous, 2003] when modeling composites with a high phase contrast using the methodology 
proposed by [Dvorak and Srinivas, 1999]. 

5.5. Influence of microstructure size. Eventually, we investigate the influence of the microstructure 
size. Figure 13 summarizes the obtained results for a moderate phase contrast and the optimal setting of 
the H-S method identified in the previous sections. A similar conclusion can be reached for the both case 
studies: for all three IjL values, the H-S method is capable of reaching the accuracy of MC confidence 
intervals for the cell length h\ approximately equal to a half of the microscopic lengthscale t. In other 
words, keeping the same number of DOFs as used to discretize the polarization problem, the accuracy 
of the method increases with the increasing l/L ratio, which is exactly an opposite trend to that of the 
classical deterministic homogenization. 




Figure 13. Influence of microstructure size; (a) statically determinate and indetermi- 
nate structures; E 2 /Ei = 5, u = 0.2, £/L = 0.05, Pq/GL 2 discretization. 
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6. Conclusions 

In the present work, the predictive capacities of numerical methods based on the Hashin-Shtrikman- 
Willis variational principles, when applied to a specific model of functionally graded materials, have been 
systematically assessed. By restricting attention to the one-dimensional setting, an extensive parametric 
study has been executed and the results of numerical schemes have been verified against reliable large- 
scale Monte-Carlo simulations. On the basis of obtained data, we are justified to state that: 

• The Hashin-Shtrikman based numerical method, when set up properly, is capable of delivering 
results with the accuracy comparable to detailed Monte Carlo simulations and, consequently, of 
outperforming the local averaging schemes. 

• When applying the Finite Element method to the solution of reference problem, the employed 
discretization has to be compatible with the numerics used to solve the polarization problem. If 
this condition is satisfied, the additional FEM-induced errors quickly become irrelevant. 

• For the discretization of the reference problem, it appears to be advantageous to combine low 
order (discontinuous) approximation of the polarization stresses with higher order quadrature 
scheme to concisely capture the heterogeneity distribution. 

• The correct choice of the reference medium has the potential to substantially decrease the error. 
Unfortunately, apart from [Dvorak and Srinivas, 1999], we fail to give any a-priory estimates of 
the optimal value for statistically non-homogeneous structures. 

• For accurate results, the characteristic cell size should be around 2-5 times smaller than the 
typical dimensions of the constituents. 

The bottleneck of the current implementation is the solution of system (32), since it leads to a fully 
populated system matrix. Fortunately, as illustrated by Figure 14, the conditioning of the polarization 
problem seems to be dominated by the phase contrast rather than the discretization of the reference 
problem, which opens the way to efficient iterative techniques. 




0.1 , 1 0.1 1 10 

(a) h/t (b) 

FIGURE 14. Sensitivity of conditioning of system matrix of the polarization problem; 
(a) statically determinate and (b) indeterminate structures; t/L = 0.05, u = 0.2, 
Po / GL2 scheme, the condition number is estimated using [Higham and Tisseur, 2000] 
algorithm. 

The next extension of the method would involve the generalization to the multi-dimensional setting. 
For the FEM-based treatment, the key aspect remains a more rigorous analysis of the combined effect of 
discretized T° operator, basis functions and integration scheme employed for the polarization problem. 
The multi-dimensional BEM approach, on the other hand, requires a careful treatment of singularities 
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of the Green function-related quantities, cf. [Prochazka and Sejnoha, 2003], which are suppressed in the 
current one-dimensional setting. Such work will be reported separately in our future publications. 
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Appendix A. Two-point probability functions 

The remaining two-point probability functions can be easily related to S\\ by exploiting the iden- 
tity (7). In particular, we obtain 

Si2(x,y) = - Su(x,y), (79) 

S21 (x, y) = Si(y) - S u (x,y), (80) 
S 2 2(x,y) = l-S 1 (x)-S 1 (y) + S u (x,y). (81) 

Appendix B. Overview of the simulation procedure 

A crude Monte-Carlo method is employed to estimate the statistics of the local fields. In particular, 
given a number of simulations N a , sampling points = yo < y\ < . . . < yjv s = L and an upper bound 
on the intensity p* > sup^rQ n p{x), the following steps are repeated for a = 1, 2, ... , N a : 

Microstructure generation: Construction of a microstructural samples is based on a two-step pro- 
cedure proposed for general Poisson processes in [Stoyan et al., 1987, Section 2.6]. First, the 
number of reference points N*(a) is determined by simulating a Poisson random variable with 
the mean p*L. The coordinates of the reference points zl(a), ^(a), . . . , Zjv*(a)( a ) tnen f° now 
from a realization of N*(a) independent random variables uniformly distributed on a closed in- 
terval [0, L]. Second, each point in the set is deleted with aprobability 1 — p yz*{a)j / p* , leading 
to a (relabeled) sequence of N p (a) particle centers z p (a). 

Solution of the one-dimensional problem: With the microstructure realization fixed, the displace- 
ment of sampling points is computed by the recursion 

( ) ( \M F m ^ - 5- m ^ rl ™ 

!iMC {y s ; oc) = umc U/s-i; a) + / — u dx, (82) 

Jy.-i E {x; a) 

where the Young modulus is provided by Equation (13) with the characteristic function xi de- 
fined as 

£ 

Xi(x;a) = 1 44> min \x — z p (a)\ > - 

p=l,2,...,N p (a) 2 

and the boundary data u(0; a) and t(0; a) determined from a generalization of the system of 
boundary equations (52)-(53). 
After completing the sampling phase, the first and second-order local statistics are assessed using the 
unbiased values 

(u)uc(y s ) ~ j-^^2u MC (y s ;a), (rlic(ys) ~ ir-—r^2{(uMc)(ys) - u M c(y s ;ac)) 

a=l a=l 

to arrive at the 7-confidence interval estimates, cf. [Rektorys, 1994, Section 34.8]: 

(u)(y s ) e [(u-)uc(ys),( u +)Mc(y s )} 

, (83) 



WMCU/sJ - t(l+-y)/2,N a -l rjrj- 1 WmCKVs) + I(l+ 7 )/2,Af Q -l 



where tp, n denotes the inverse of the Student t distribution function for value j3 and n degrees of freedom. 
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The reference results reported in Section 5 correspond to the values obtained for N a = 100, 000 
simulations, the confidence level 7 = 99.9%, 101 equidistant sampling points and the integral (82) 
evaluated with an adaptive Simpson quadrature [Gander and Gautschi, 2000] with the relative tolerance 
set to 10" 6 . 
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